Dynamical density-density correlations in the one-dimensional Bose gas 
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The zero-temperature dynamical structure factor of the one- dimensional Bose gas with delta- 
function interaction (Lieb-Liniger model) is computed as a function of momentum and frequency 
using a hybrid theoretical/numerical method based on the exact Bethe Ansatz solution. This allows 
to interpolate continuously between the weakly-coupled Thomas-Fermi and strongly-coupled Tonks- 
Girardeau regimes. The results should be experimentally accessible with Bragg spectroscopy. 



The physics of low dimensional atomic systems 
presents very special features as compared to the three- 
dimensional case. As the temperature is lowered, a uni- 
form gas of bosons in three dimensions will undergo 
a transition to a Bose-Einstein condensate (BEC)^; in 
the one-dimensional case low-energy fluctuations pre- 
vent long-range order. For trapped gases, the situation 
changes and three regimes become possible in ID^: true 
condensate, quasicondensate, and a strongly-interacting 
regime, with BEC limited to extremely small interaction 
between particles. Trapped ID gases are now accessible 
experimentally^i^ in all regimes, the most challenging to 
obtain being the strongly-interacting case^ii, which can 
survive without fast decay due to a reduced three-body 
recombination rate^ (a consequence of fermionization) . 

A natural starting point for the theoretical descrip- 
tion of one-dimensional atomic gases in this last regime 
is provided by bosons with delta-function interaction (the 
Lieb-Liniger modelAS) , whose Hamiltonian is given by 
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in which c > is the coupling constant, and the sum is 
over pairs (we have put h = 1 = 2m for simplicity). For 
definiteness, we consider a system of length L with peri- 
odic boundary conditions. In the thermodynamic limit, 
the physics of the model depends on a single parame- 
ter 7 = c/n where n = N/L is the particle density. In 
ID, in stark contrast to higher dimensions, low densities 
lead one to the strong-coupling regime of impenetrable 
bosons, know as the Tonks-Girardea u^^'^^ limit. 

Although equilibrium thermodynamic properties of the 
Lieb-Liniger model are accessible via the Bethe Ansatz^, 
dynamical objects such as correlation functions cannot be 
readily obtained with this scheme. For example, the zero- 
temperature density-density correlation function (writ- 
ten here in Fourier space, where it is also known as the 
dynamical structure factor (DSF)), 



S{k,uj) = 



dx I dt e 



{p{x,t)p{0,0)) (2) 



(in which p{x) — ^(■^ ~ -^i)) '^o^ resisted 

all efforts towards an exact computation. The present 
paper presents a reliable and efficient method for com- 
puting this, based on mixing integrability and numerics. 



Many approximate theoretical schemes have been de- 
veloped to tackle this issue. In the BEC regime, Bogoli- 
ubov theory can be used in conjunction with local den- 
sity, impulse or eikonal approximations^^. Specifically 
in ID, an effective harmonic fiuid approach (Luttinger 
liquid theory)^^ can be used to obtain information on 
the asymptotics of static and dynamical correlation func- 
tions at zero and nonzero temperatur o^^'^^i^^ . Inversely 
to asymptotics, a small distance Taylor expansion was 
also proposedi^. Yet another possibility is to exploit an 
exact fermion mapping, and use the Hartree-Fock and 
generalized random phase approximation to get dynam- 
ical correlators near the Tonks-Girardeau limit in a I/7 
expansior^. Quantum Monte Carlo has been used to 
study this limiii^i^ and to numerically obtain the pair 
distribution function and static structure factor—. How- 
ever, up to now, there is no overall reliable method for 
obtaining the full momentum and frequency dependence 
of the DSF. 

In view of the integrability of ([T]), one could expect 
to obtain nonperturbative results for objects such as ^ . 
Much recent progress on the computation of correlation 
functions for this and other ID quantum integrable mod- 
els has in fact been achieved through the Algebraic Bethe 
Ansata^^i^i^. In this paper, we wish to present a novel 
method for obtaining dynamical correlation functions of 
model ([!]), which is based on these developments. We 
will obtain the dynamical structure factor for finite but 
large systems, starting directly from the Bethe Ansatz 
solution, in a way which is reminiscent of recent work by 
one of us on dynamical spin-spin correlation functions in 
Heisenberg magnets^^. In particular, the momentum and 
frequency dependence of the DSF is fully characterized 
by our approach. All our results are presented in Figures 
1-3. The static structure factor S{k) = / |^S'(fc,w) is 
also obtained as a subset of our results. The DSF itself 
is experimentally accessible through Fourier sampling of 
time-of-flight images^! or through Bragg spectroscopy^^. 

By inserting a summation over intermediate states, 
S{k,uj) is transformed into a smu of matrix elements of 
the density operator in the basis of Bethe eigenstates \a), 



where pu = Z^jLi e 
term of this sum are fully determined by the Algebraic 



||pfe|a)p5(w-^„+^o) (3) 



All of the elements in each 
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FIG. 1: Density plots of the dynamical structure factor as 
a function of interaction. The horizontal axis is momentum 
(running up to Akp), and the vertical axis represents energy 
transfer. Data obtained from systems of length L — 100 with 
N = 100 and 7 = 0.25, 1, 5, 10, 20 and 100. 

Bethe Ansatz, for a finite system with specified boundary 
conditions (we consider periodic systems in the present 
work). The Fock space is spanned by the set of Bethe 
wavef unctions, each fully determined by a set of N ra- 
pidities {Xj}, solution to the Bethe equations 

A,- + - > 2 arctan = — (4) 

fc c L 

in which Ij are half-odd integers (integers) for even (odd) 
iV. All solutions to these Bethe equations are real. Each 
set of distinct quantum numbers {Ij}, Ij ^ Ik '^i j ^ k 
defines a Bethe eigenstate participating in the sum ([3]). 
The energy and momentum of such a state are given by 
E = and k = Y.jXj = xEj^- The ground- 

state itself is obtained from the set {I^}, with 7° = -^^j-^- — 
j, j = 1, TV. The wavefunction of an eigenstate is given 
by the Bethe Ansatz, and its norm by the determinant 
of the Gaudin matri x^^i^" . 

Matrix elements of the density operator in the basis 
of Bethe eigenstates were calculated with the Algebraic 
Bethe Ansatz in [2^ . They are given by the determinant 
of a matrix whose entries are rational functions of the 
rapidities of the two eigenstates involved. For the sake 
of brevity we do not reproduce these expressions here. 

What remains to be performed is the actual summation 
over intermediate states in ([3]). From this step onwards. 



TABLE L /-sum rule saturation percentage achieved at the 
two representative values of momentum used in Fig. 2. The 
method converges fastest in the strongly-interacting regime, 
the most difficult regime being intermediate interactions. 



7 


0.25 


1 


5 


10 


20 


100 


kp 


99.52 


99.44 


99.48 


99.45 


99.80 


99.97 


2kF 


99.34 


99.11 


97.49 


98.21 


99.35 


99.90 



TABLE II: Sound velocities for the interaction values consid- 
ered here. A comparison is made between the finite-size value, 
and the infinite-size one. 
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everything is done numerically. The Fock space of inter- 
mediate states is scanned by navigating through choices 
of sets of quantum numbers. For each individual inter- 
mediate state, the Bethe equations are solved, and the 
matrix element is computed. To obtain smooth curves in 
energy, the energy delta function in ([3]) is broadened to a 
width equal to a multiple of the typical energy level spac- 
ing. The contribution to the dynamical structure factor 
sum is tallied until good convergence has been achieved. 
This is quantified by evaluating the /-sum rule, 

j'^.s{k,.)^^e. (5) 

Since this is skewed towards high energy, and in view of 
the ordering of states in the scanning we perform (typi- 
cally going from low-energy intermediate states to higher- 
energy ones), the saturation level of this sum rule repre- 
sents a lower bound for the saturation of S{k,Lij) itself. 

It is useful to recall here the nature of excitations 
in the Lieb-Liniger model, which come in two typesi^, 
Type I ("particles") and Type II ( "holes" )3i. Type I 
are Bogoliubov-like quasiparticles that exist for any mo- 
mentum, and represent states with one quantum number 
displaced outside the ground-state interval. Their dis- 
persion relation is described in the thermodynamic limit 
by an integral equation, yielding a curve contained be- 
tween the asymptotic limits e/(fc) = fc^ at 7 = and 
ti{k) = k'^ + 2TTn\k\ for 7 — > 00. Type II excitations 
are holes in the ground-state distribution, and do not 
appear in Bogoliubov theory. They exist in the interval 
fc| < kp = 7rn, and their dispersion relation coincides 
with the lower threshold of the DSF. Low-energy Umk- 
lapp modes aX k = 2kp can be understood as excitations 
going from one side of the Fermi surface to the other. 
Both Type I and Type II dispersion relations approach 
fc — > with a slope equal to the velocity of sound Vs{'^)- 

A remarkable feature in this method comes from the 
fact that it is not necessary to scan through the whole 
Fock space to get good saturation. Contributions from 
intermediate states with up to only a handful of particles 



FIG. 2: Fixed momentum plots of the dynamical structure factor for five representative values of the interaction parameter 7 
obtained using L = 100 and A'^ — 100 (continuous curves) and L = 80 and A'' = 80 (dots), such that kp = n in all cases. The 
energy (5-function in equation ((3]) is a Gaussian of width w = 0.3. 
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FIG. 3: Static structure factor for four representative values 
of the interaction parameter 7. The three values of system 
size illustrate the smallness of finite-size effects. All curves 
saturate towards the limit S{k —> 00) = 1. The dashed lines 
are the small momentum asymptotic predictions S{k) ~ k/vs- 
Deviations in the 7 = 5 curve near 2fe_F are commensurate 
with lack of saturation of the f-sum rule there (see Table I). 



are sufficient to achieve extremely good accuracy. The 
saturation of the /-sum rule ([5]) achieved for the TV = 100 
curves in Fig. 2 is summarized in Table 1 (for lower values 
of N, even better saturation is achieved). Our algorithm 
is designed to recursively hunt for the most important 
terms in the multiparticle sum, in decreasing order of 
contribution to the DSF, and therefore to maximize the 
efficiency of the available numerical resources. Higher ac- 
curacy is obtained by using more computational time to 
include the contributions from more intermediate states. 

Fig. 1 shows density plots for the DSF obtained with 
our method, for different values of 7. We used unit den- 
sity n = 1 with L = 100 and = 100 (this compares 



with experimental particle numbers^i). For small 7, the 
DSF is essentially a single delta peak centered on the 
Type I dispersion relation, S{k,uj) = T7T(k)^^'^ ~ ^^(^))- 
As 7 increases, the peak flattens but remains near the 
upper two-particle boundary. In particular, this con- 
firms that the peak near the lower boundary observed 
in first-order RPA is an artefact of that method^S. In- 
creasing 7 further towards the Tonks-Girardeau limit, 
the DSF approaches a constant over a finite frequency 
interval for a given momentum. All of this is illustrated 
more specifically for two specific values of momentum in 
Fig. 2. Fig. 3 shows S'(fc), which qualitatively fits with 
quantum Monte Carlo results2^. At small momentum, 
we recover the prediction S{\z) ~ fc/wg (see for exam- 
ple [l^l). The curve for the largest interaction value, 
7 — 100, also clearly approaches the well-known result 
lim-y^oo 5'(fc) = kjlkp for fc < 2kp. Moreover, Fig. 1 
shows that low-energy contributions near 2fci?, which rep- 
resent superfluidity-breaking Umklapp modes, are only 
important for large 7. 

For all values of 7, the signal mostly lies between 
the two-particle continuum defined by convolution of the 
Type I and Type II dispersion relations. For general fc, 
all the signal in fact lies strictly above the Type II dis- 
persion relation modulo 2kp translations, since there are 
no lower-energy states available. The upper two-particle 
bound is however not robust: multiparticle contributions 
give a nonzero signal at (in principle arbitrarily) large en- 
ergy (coming from e.g. states with two or more particles 
of large but opposite momentum). In practice, however, 
the data shows that the onset of the DSF at finite 7 is 
followed by a sharp peak around the Type I dispersion, 
followed by a rapid decrease. We believe that in the ther- 
modynamic limit, contributions from intermediate states 
with higher particle numbers smoothen the upper thresh- 
old into a high-frequency tail, as is the case for the corre- 
sponding correlators in quantum spin chains'^'*. The only 
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exception is the Tonks-Girardeau limit, where both the 
lower and upper thresholds remain sharp. 

To quantify finite-size effects in our results, we have 
included data for = 80 in Fig. 2 and for N = 60, 80 in 
Fig. 3 and given values for the effective ^8(7) (the con- 
formal exponent 77, such that S{2kF,i^) — oj^^'^ at small 
Lu, is given by 77 = Ann/vs) obtained for iV = 100 com- 
pared to the thermodynamic one in Table [TTl Theoretical 
considerations based on Bethe Ansatz expressions for cor- 
relation functions^, similar to those used here, predict 
finite-size corrections of order 1/N. The precise form 
of these 1/N corrections then depends on the specific 
boundary conditions used, but we believe that our results 
are close enough to the thermodynamic limit to make 
the choice of boundary conditions immaterial. Specifi- 
cally, in Fig. 2, the only observable effect of increasing 
system size is a very slight shift of the main peaks of the 
DSF towards higher energy (for clarity we do not plot the 

= 80 results in the right-hand figure, but they show 
the same behaviour as in the left-hand one). The static 
structure factor plotted in Fig. 3 shows essentially no 
change for the different values of N given. The variations 
observed fit comfortably within the deviation from per- 
fection of the sum rule saturation achieved, which is the 
actual determining factor in the quality of our results. 
More important for theory is the question of a confin- 



ing potential2^i^i^, which is present in experiments but 
breaks the integrability of the Lieb-Liniger model. We 
expect that experiments on large enough systems would 
on the other hand show correlations approaching those of 
a pure Lieb-Liniger model similar to those obtained here, 
or that experiments could be done in box-like geometries, 
where a variation of our method would be applicable. 

Summarizing, we have computed the frequency- and 
momentum-dependent dynamical density-density corre- 
lation function of the one-dimensional interacting Bose 
gas (Lieb-Liniger model) for systems with finite num- 
bers of particles, using a Bethe Ansatz-based numeri- 
cal method. This goes beyond other available methods 
in offering a full characterization of the momentum and 
frequency dependence of the dynamical structure factor, 
provides a firm testing standard for other methods, and 
opens the way to many possible extensions (other correla- 
tors, systems with mixed statistics, finite temperatures) 
on which we will report in future publications. 
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